########Main Analysis Replication for Trussler BJPS 2019#########
rm(list=ls())
library(sandwich)
library(foreign)
library(stargazer)
library(car)
library(plm)
library (sp)
library (rgdal)
library (raster)
library(maptools)
library(RColorBrewer)

#Change to working directory in which you have the data files
setwd("C:/Google Drive/Information and Accountability/Incumbency Bias Nationalization/Trussler BJPS Replication/ReplicationData")


#####Incumbency Advantage Over Time#####
jacob <- read.dta("JOPrepfile1.dta")


#Calculate incumbency advantage by year (Jacobson Replication)

year <- unique(jacob$year)
year <- year[year>=1960 & year<=2014]
year <- year[!is.na(year)]
incumbency.advantage.south <- NA
ia.se.south <- NA
incumbency.advantage <- NA
ia.se <- NA


for(i in 1:length(year)){
  m <-   lm(dv ~ ptynow + dpres +inc3, data=jacob[jacob$year==year[i] & jacob$south==1,])
  incumbency.advantage.south[i] <- coef(m)["inc3"]
  ia.se.south[i] <- sqrt(vcov(m)["inc3","inc3"])
  m <-   lm(dv ~ ptynow + dpres +inc3, data=jacob[jacob$year==year[i]& jacob$south==0,])
  incumbency.advantage[i] <- coef(m)["inc3"]
  ia.se[i] <- sqrt(vcov(m)["inc3","inc3"])
}

incumbency.advantage

#Create Figure One: Incumbency Advantage over Time for South and Non-South
pdf(file="./MainFigures/Figure1.pdf", width = 11, height = 8)
plot(year, incumbency.advantage, type="b", pch=16, lwd=2, ylim=c(-0,16), xlim=c(1960,2020),
     axes=F, xlab="", ylab="Incumbency Advantage")
points(year, incumbency.advantage.south , pch=17, type="b", col="gray60", lwd=2)
axis(side=1, at=seq(1960,2020,10))
axis(side=2, at=seq(0,15,5))
legend("topright", c("Non-South","South"), pch=c(16,17), col=c("black","gray60"))
dev.off()





######Broadband Descriptives######

#Load broadband data
load(file="Broadband by Congressional District.Rdata")



#Make Broadband Map for 2002 
cd.bb$cdfip <- substr(cd.bb$unique,5,8)


cd02 <- cd.bb$cdfip[cd.bb$year==2002]
providers02 <- cd.bb$providers[cd.bb$year==2002]
prov.data02 <- cbind.data.frame(cd02, providers02)




districts.108 <- readOGR(dsn = "cd99_108_shp", layer="cd99_108")

head(districts.108)
table(districts.108$STATE)
table(districts.108$CD)

districts.108$CD[districts.108$CD=="00"] <- "01"

districts.108$cdfip <- paste(districts.108$STATE,districts.108$CD,sep="")

districts.108 <- merge(districts.108,prov.data02, by.x="cdfip",by.y="cd02")


names(districts.108)

districts.108 <- districts.108[!is.na(districts.108$STATE),]
districts.108 <- districts.108[districts.108$STATE!=72,]
districts.108 <- districts.108[districts.108$STATE!=15,]
districts.108 <- districts.108[districts.108$STATE!="02",]

my.palette <- brewer.pal(n = 6, name = "Greys")
my.palette <- my.palette[-5]



#National Map

c(2,4,6,8,10,12)

#Just for this map collapse top of scale to 12
districts.108$providers02[districts.108$providers02>=12] <- 11.99

##Uncomment to make map (computer intensive)

#pdf(file="./MainFigures/NationalBroadbandMaplegend.pdf", height=10, width=10)
#spplot(districts.108, "providers02",col.regions = my.palette,at=c(2,4,6,8,10,12), col ="gray60")
#dev.off()



#Levels and change in Broadband Providers

cd.bb$cdfip <- as.character(cd.bb$cdfip)
cd.bb <- cd.bb[cd.bb$year>=2002,]

cds8 <- unique(cd.bb$cdfip[cd.bb$year==2008]) 
cds2 <- unique(cd.bb$cdfip[cd.bb$year==2002]) 
cds <- cds8[cds8 %in% cds2]

cd.bb$prov.change <- NA
for(i in 1:length(cds)){
  cd.bb$prov.change[cd.bb$cdfip==cds[i] & cd.bb$year==2004] <- cd.bb$providers[cd.bb$cdfip==cds[i] & cd.bb$year==2004] - cd.bb$providers[cd.bb$cdfip==cds[i] & cd.bb$year==2002]
  
  cd.bb$prov.change[cd.bb$cdfip==cds[i] & cd.bb$year==2006] <- cd.bb$providers[cd.bb$cdfip==cds[i] & cd.bb$year==2006] - cd.bb$providers[cd.bb$cdfip==cds[i] & cd.bb$year==2004]
  
  cd.bb$prov.change[cd.bb$cdfip==cds[i] & cd.bb$year==2008] <- cd.bb$providers[cd.bb$cdfip==cds[i] & cd.bb$year==2008]  - cd.bb$providers[cd.bb$cdfip==cds[i] & cd.bb$year==2006] 
}

pdf(file="./MainFigures/ProvidersLevels.pdf")
boxplot(cd.bb$providers ~cd.bb$year, ylim=c(0,20), outline=F, ylab="Number of Broadband Providers", axes=F)
axis(side=1, at=seq(1,4,1), labels=c( "2002","2004","2006","2008"))
axis(side=2, at=seq(0,20,5), las=2)
dev.off()

pdf(file="./MainFigures/ProvidersChange.pdf")
boxplot(cd.bb$prov.change ~cd.bb$year, outline=F, axes=F, ylim=c(-3,8), ylab="Change in Broadband Providers")
axis(side=1, at=seq(1,3,1), labels=c( "2002-2004","2004-2006","2006-2008"))
axis(side=2, at=seq(-2,8,2), las=2)
abline(h=0, lty=2)
dev.off()

###### Main Analysis ######
load(file="AggregateNatInc.Rdata")

###Incumbency Advantage


m.incumbency1 <- plm(dv ~ factor(year) + ptynow + dpres +inc3   ,
                     data=cd, model="within", index=c("cdfip","year"))
m.incumbency1.vcov <- vcovHC(m.incumbency1, type="HC0", cluster="group", method = "arellano")
m.incumbency1.se <- sqrt(diag(m.incumbency1.vcov))
m.incumbency2 <- plm(dv ~factor(year) + ptynow + dpres + inc3*log(providers) + 
                       log(med.income) + log(pop) ,
                     data=cd, model="within", index=c("cdfip","year"))
m.incumbency2.vcov <- vcovHC(m.incumbency2, type="HC0", cluster="group", method = "arellano")
m.incumbency2.se <- sqrt(diag(m.incumbency2.vcov))
summary(m.incumbency2)


#Creation of Table 1
stargazer(m.incumbency1,m.incumbency2, omit=c("cdfip", "year"),
          se=list(m.incumbency1.se, m.incumbency2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          covariate.labels = c("Party Now$_{kt}$","% Democratic Vote for President$_{kt}$","Incumbency$_{kt}$",
                               "log(Broadband Providers$_{kt}$)","log(Median Income$_{kt}$)",
                               "log(Population$_{kt}$)", 
                               "Incumbency$_{kt}$*log(Broadband Providers$_{kt}$)"),
          dep.var.labels   = "% Democratic Vote for House$_{kt}$",
          model.numbers = F,
          add.lines = list(c("District F.E.", "Yes","Yes"),
                           c("Election Year F.E.","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incumbency")


#Partisan Effects on Incumbents Models

m.incwin.pres1 <- plm(inc.vote ~ factor(year)*factor(party) +   pres.vote +
                        log(pop)  + log(med.income)  , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres1.vcov <-  vcovHC(m.incwin.pres1, type="HC0", cluster="group", method = "arellano")
m.incwin.pres1.se <- sqrt(diag(m.incwin.pres1.vcov))
summary(m.incwin.pres1)


m.incwin.pres2 <- plm(inc.vote ~ factor(year)*factor(party)  +  log(providers)*pres.vote +
                        log(pop)  + log(med.income)  , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.pres2.vcov <-  vcovHC(m.incwin.pres2, type="HC0", cluster="group", method = "arellano")
m.incwin.pres2.se <- sqrt(diag(m.incwin.pres2.vcov))
summary(m.incwin.pres2)



m.incwin.party1 <- plm(inc.vote ~ factor(year)*factor(party)    + 
                         log(pop)  + log(med.income)  , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party1.vcov <-  vcovHC(m.incwin.party1, type="HC0", cluster="group", method = "arellano")
m.incwin.party1.se <- sqrt(diag(m.incwin.party1.vcov))


m.incwin.party2 <- plm(inc.vote ~factor(year)*factor(party)   +  log(providers)*party.unity  + 
                         log(pop)  + log(med.income) , data=cd, model="within", index=c("icpsr.id","year"))
m.incwin.party2.vcov <-  vcovHC(m.incwin.party2, type="HC0", cluster="group", method = "arellano")
m.incwin.party2.se <- sqrt(diag(m.incwin.party2.vcov))
summary(m.incwin.party2)


save(m.incumbency2, m.incumbency2.vcov,
     m.incwin.pres2, m.incwin.pres2.vcov,
     m.incwin.party2, m.incwin.party2.vcov,
     file="MainModels.Rdata")

stargazer(m.incwin.pres1, m.incwin.pres2,m.incwin.party1, m.incwin.party2, omit=c("icpsr.id","name","year","ptynow"),
          order=c("pres.vote", "party.unity","log(providers)", "log(providers):pres.vote","log(providers):party.unity"),
          se=list(m.incwin.pres1.se, m.incwin.pres2.se,m.incwin.party1.se, m.incwin.party2.se), type="text", report='vc*p')

#Creation of Table 2
stargazer(m.incwin.pres1, m.incwin.pres2,m.incwin.party1, m.incwin.party2, omit=c("icpsr.id","name","year","ptynow"),
          se=list(m.incwin.pres1.se, m.incwin.pres2.se,m.incwin.party1.se, m.incwin.party2.se),star.char = c("*"),
          star.cutoffs = c(0.05), style="ajps",omit.stat = c("adj.rsq", "f"), digits=2,
          order=c("pres.vote", "party.unity","log(providers)", "log(providers):pres.vote",
                  "log(providers):party.unity"),
          covariate.labels = c("Same Party Presidential Vote$_{jt}$",
                               "Same Party Presidential Vote$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "Party Unity$_{jt}$",
                               "Party Unity$_{jt}$*log(Broadband Providers$_{jt}$)",
                               "log(Broadband Providers$_{jt}$)", 
                               "log(Population$_{jt}$)",
                               "log(Median Income$_{jt}$)" ),
          dep.var.labels   = "% Incumbent Vote$_{jt}$",
          model.numbers = F,
          add.lines = list(c("Member F.E.", "Yes","Yes","Yes","Yes"),
                           c("Election Year*Party F.E.","Yes","Yes","Yes","Yes")),
          notes        = "$^{*}$p $<$ .05; Cluster-Robust Standard Errors", 
          notes.append = FALSE,
          title= "Effect of Communication Environment on Incumbency Advantage",
          label="incvote")




providers.sample  <- sample(cd$providers,500)
Rmarginsln <- function(model,xname,xzname,vcov.matrix,levels, ci=.95, latex=F){
  # Generate Marginal Effects
  m_effect <- coef(model)[xname]+coef(model)[xzname]*levels
  #Variances and Covariances
  v_x <- vcov.matrix[xname,xname]
  v_xz <- vcov.matrix[xzname,xzname]
  cov <- vcov.matrix[xname,xzname]
  #Standard Errors
  se <- sqrt(v_x + (levels^2)*v_xz+2*levels*cov)
  #T value for 95%CI
  if (ci==.95){
    t <- qt(0.025,model$df)
  }
  #T value for 90%CI
  if (ci==.9){
    t <- qt(0.05,model$df)
  }
  #Confidence Bounds
  m_upper <- m_effect + se*t
  m_lower <- m_effect - se*t
  # Remove Flotsom and Jetson
  #Printing Table
  table <- cbind(levels,m_effect,se,m_upper,m_lower)
  if (latex==T){
    library(xtable)
    print(xtable(table),include.rownames=FALSE)
  }
  if (latex==F){
    return(table)
  }
}

lnproviders <- log(seq(.1,20,.01))
margins <- Rmarginsln(m.incwin.pres2,"pres.vote", "log(providers):pres.vote", m.incwin.pres2.vcov, levels = lnproviders)

Rmarginsln(m.incwin.pres2,"pres.vote", "log(providers):pres.vote", m.incwin.pres2.vcov, levels = log(median(cd$providers,na.rm=T)))

pdf("./MainFigures/PresVote.pdf", height=6, width=8)
plot(exp(margins[,1]), margins[,2], type="l", ylim=c(-.6,.6), xlim=c(2,16), axes=F,
     xlab="Broadband Providers", ylab="Effect of District Presidential Vote on Incumbent Margin", lwd=2)
lines(exp(margins[,1]), margins[,4], lty=2)
lines(exp(margins[,1]), margins[,5], lty=2)
rug(providers.sample)
abline(h=0, lty=2, col="gray60")
axis(side=1, at=seq(2,16,2))
axis(side=2, at=c(-.6,-.4,-.2,0,.2,.4,.6), las=2)
dev.off()


lnproviders <- log(seq(.1,20,.01))
Rmarginsln(m.incwin.party2,"party.unity", "log(providers):party.unity", m.incwin.party2.vcov, levels = log(median(cd$providers,na.rm=T)))

margins <- Rmarginsln(m.incwin.party2,"party.unity", "log(providers):party.unity", m.incwin.party2.vcov, levels = lnproviders)

pdf("./MainFigures/PartyUnity.pdf", height=6, width=8)
plot(exp(margins[,1]), margins[,2], type="l", ylim=c(-.6,.6), xlim=c(2,16), axes=F,
     xlab="Broadband Providers", ylab="Effect of Party Unity on Incumbent Margin", lwd=2)
lines(exp(margins[,1]), margins[,4], lty=2)
lines(exp(margins[,1]), margins[,5], lty=2)
rug(providers.sample)
abline(h=0, lty=2, col="gray60")
axis(side=1, at=seq(2,16,2))
axis(side=2, at=seq(-6,.6,.2), las=2)
dev.off()


